{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfafc213-e317-4e1b-a977-9d91c188a29d",
   "metadata": {},
   "outputs": [],
   "source": [
    "suppressMessages(library(ArchR))\n",
    "suppressMessages(library(Seurat))\n",
    "suppressMessages(library(Signac))\n",
    "suppressMessages(library(cowplot))\n",
    "suppressMessages(library(dplyr))\n",
    "suppressMessages(library(tidyr))\n",
    "suppressMessages(library(ComplexHeatmap))\n",
    "suppressMessages(library(RColorBrewer))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0694d447-00f5-4767-8de5-ea66de672e2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj.atac <- readRDS(\"../data/VisiumHeart/snATAC.annotated.Rds\")\n",
    "obj.atac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6c31ca5-4d0d-43b4-ae84-b417be58048f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## add gene activity matrix\n",
    "gene.activities <- readRDS(\"../data/VisiumHeart/GeneScoreMatrix.Rds\")\n",
    "\n",
    "gene.activities <- gene.activities[, colnames(obj.atac)]\n",
    "\n",
    "obj.atac[['GeneActivity']] <- CreateAssayObject(counts = gene.activities)\n",
    "\n",
    "DefaultAssay(obj.atac) <- \"GeneActivity\"\n",
    "\n",
    "obj.atac <- obj.atac %>% \n",
    "        NormalizeData() %>%\n",
    "        FindVariableFeatures() %>%\n",
    "        ScaleData() %>%\n",
    "        RunPCA(verbose = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "079c4590-b13d-4edb-88d3-9e699b4e1b68",
   "metadata": {},
   "outputs": [],
   "source": [
    "Idents(obj.atac) <- \"cell_type\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca0f4556-bed3-459c-aad6-0f202b3f19ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "all.markers <- FindAllMarkers(obj.atac, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "285a521a-77bc-4285-ad21-3670ccb56cd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "all.markers$cluster <- factor(all.markers$cluster, levels = c(\"CM\", \"Fib\", \"Endo\", \"Lymphoid\", \"Myeloid\", \"Pericyte\", \"vSMCs\", \"Neuronal\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a76d80cd-94dd-4f26-a4b8-cb6dc0cc0e2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "all.markers <- all.markers[order(all.markers$cluster), ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fc19293-0bc3-4e2e-b6cc-f38015db4064",
   "metadata": {},
   "outputs": [],
   "source": [
    "df <- all.markers %>%\n",
    "    group_by(cluster) %>%\n",
    "    slice_max(n = 10, order_by = avg_log2FC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db8a5caa-1223-4df3-a421-7002da63e11b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## get average gene activity score\n",
    "avg <- AverageExpression(\n",
    "  obj.atac,\n",
    "  assays = \"GeneActivity\",\n",
    "  features = df$gene,\n",
    "  return.seurat = FALSE,\n",
    "  group.by = \"cell_type\",\n",
    "  slot = \"data\",\n",
    "  verbose = TRUE,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "627ed0f2-ec16-4007-84a9-2fa0217fba55",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot <- avg$GeneActivity\n",
    "\n",
    "head(df_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "803e1e54-fb1e-4095-aa62-08a60bc82ebe",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot_scale <- t(apply(df_plot, 1, scale))\n",
    "colnames(df_plot_scale) <- colnames(df_plot)\n",
    "rownames(df_plot_scale) <- rownames(df_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79c4dd13-8824-40ab-8f46-d08abc8180bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(df_plot_scale)\n",
    "\n",
    "df_plot_scale <- df_plot_scale[, c(\"CM\", \"Fib\", \"Endo\", \"Lymphoid\", \"Myeloid\", \"Pericyte\", \"vSMCs\", \"Neuronal\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8244014d-4aee-4ea7-a150-ff2408fdf7c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot_scale <- df_plot_scale[rev(rownames(df_plot_scale)), ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b3cc272-55b2-4412-a33f-41844c596445",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 4, repr.plot.height = 14)\n",
    "\n",
    "p <- Heatmap(as.matrix(df_plot_scale),\n",
    "             name = \"Gene Activity\",\n",
    "             cluster_columns = FALSE,\n",
    "             cluster_rows = FALSE,\n",
    "             show_row_names = TRUE,\n",
    "             rect_gp = gpar(col = \"black\", lwd = 0.5),\n",
    "             col = ArchR::paletteContinuous()\n",
    "            )\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "142bf51e-67ce-4ba7-bc61-2f98c6b08ee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "saveRDS(df_plot_scale, \"../data/VisiumHeart/marker_heatmap.Rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf49fb87-ced1-4331-99ce-f82f7e6db7d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Session information\n",
    "sessionInfo()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
